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Study of Liapunov Exponents and the Reversibility of Molecular 
Dynamics Algorithms 
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We study the question of lack of reversibility and the chaotic nature of the equations of motion in numerical 
simulations of lattice QCD. 
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1. Introduction 

Molecular dynamics algorithms like the Hybrid 
Monte Carlo and the Kramers equation 
algorithm have been playing a major role in nu- 
merical simulations of QCD on a Euclidean lat- 
tice. In order for these algorithms to fulfill the 
detailed balance condition, the classical motion 
governed by the set of Hamilton's equations for 
the system should be reversible. However, as no- 
ticed some time ago this reversibility condi- 
tion is violated due to the round-off errors in the 
numerical integration of the equations of motion. 
Recently, it was pointed out |^ that, due to the 
chaotic nature of the equations of motion, the 
round-off errors in the integration which violate 
the reversibility condition get magnified exponen- 
tially with a positive Liapunov exponent ly. 

We will focus our discussion on the equations 
that arise in the simulations of Wilson QCD with 
two flavors of quarks with degenerate masses. Pe- 
riodic boundary conditions have been taken for all 
the fields. The full partition function for Wilson 
QCD is given by, 



[-Sg-cj,^Q-^[U]cb) ,(1) 



with Sg being the Wilson plaquette action. The 
fermion matrix (5[C/] is a hermitian sparse matrix. 

In molecular dynamics algorithms , one in- 
troduces the Hamiltonian: 
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-I- SeffiUx,f_„ 0^ 



(2) 



where H^^fi is the Gaussian distributed momen- 
tum conjugate to the gauge field C/x,p and takes 
values in su(3), the Lie algebra of SU{3). We 
have also introduced the effective action Seff = 
Sg + 4>^Q~^(f>. Then the gauge fields and its cor- 
responding momenta are updated according to 
Hamilton's equations of motion: 



iHx.fj. = [Ux.fj.Fx^^i]T.A. , (3) 



where the symbol [• • ■]t.a. stands for taking the 
traceless antihermitian part of the matrix |^ and 
the quantity Ux,^iFx,^i is the total force associated 
with the link Ux,ti- The dot on a field variable in- 
dicates the derivative with respect to Monte Carlo 
time. Eq. (||) defines a Hamilton flow in a phase 
space manifold which is a direct product of AL^T 
factors of SU{'i) and su(3). 



2. Liapunov Exponents 

For any flow described by a set of first-order 
autonomous differential equations, the concept of 
Liapunov exponents [^,0 could be introduced lo- 
cally at each point in the phase space manifold. 
They describe the mean exponential rate of di- 
vergence of two nearby trajectories. 

One is thus led to study the time evolution 
of tangent vectors. We take another differen- 
tial of eq. (|^) and denote dHx^^ and dXx^^ = 
~^^x,li'^Ux,tJ. as the corresponding tangent vec- 
tors, we have: 



idH, 



x,^ 



^x,fj.'^^x,i.iUx,fj, , 
= d[Ux^p,Fx^f^]T.A. ■ 



(4) 
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We introduce the norm in tangent space as: 
D\t) = ^tr(dij2^^(r) + dXl^ir)) . (5) 



The Liapunov exponent can be defined as: 

h' = lim lim — log — . (6) 

D(0)^OT^oo T ^ D{0) ^ ' 

The exponent with the largest real part is called 
the leading Liapunov exponent. 

The numerical calculation of Liapunov expo- 
nents of a given flow can be done straightfor- 
wardly [||. Starting at r = 0, we have some 
initial value of D^{Q) = Y..,^.tr{dHl^ + dXl^) 
for a given initial tangent vector {dHx,fj.,dXx^fj,). 
We then integrate one step in time with step size 
St using the leapfrog scheme. Now, we eval- 
uate the norm D{5t) of the new tangent vec- 
tor and store this information. Next, we rescale 
the new tangent vector such that its norm is 
still equal to D{0). Repeating this for N^d 
steps, we get a sequence of values for the norm: 
D{0),D{ST),---D{N^dT). It can be shown § 
that the average 
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D{k5T) 



(7) 



is approaching the leading Liapunov exponent 
when n ^ oo. Its value is independent of the 
value of 5t, as long as 5t is not too large. 

3. Liapunov Exponents for Simulations of 
QCD 

We have studied the leading Liapunov expo- 
nents for various values of (3 and k on 4'' lattices. 
In Fig. 1^ (a), we plot the exponents as a func- 
tion of (3 for the pure SU (3) gauge theory. ^ We 
see from this figure that the Liapunov exponents 
show a significant /3-dependence. It has been con- 
jectured recently 1^ that this dependence might 
be related to the correlation length of the theory. 



^Similar studies have also been done by the authors of 
ref . fel . Complete consistent values of the exponent up to 
/3 = 10 have been obtained. Note that there is a factor of 
\/2 difference between their normalization of Monte Carlo 
time and ours. 




Figure 1. The average leading Liapunov expo- 
nents for pure SU{?>) gauge theory (a) and full 
QCD (b) as a function of the gauge coupling (3. 
In (b) , exponents for different values of k are rep- 
resented by different symbols. The triangles cor- 
respond to K — 0.13 and the hexagons correspond 
to K = 0.16. 



However, this hypothesis is subject to further in- 
vestigation. 

We have also performed similar investigations 
of the Liapunov exponents for full QCD with dy- 
namical fermions in the parameter range 5.6 < 
(3 < 8.0 and 0.13 <k< 1.6. In Fig. 0(b), we show 
the final result of the exponents as a function of /3 
and K. For each value of f3, we have taken 4 differ- 
ent values of k, namely 0.13, 0.14, 0.15 and 0.16. 
Due to their weak dependence on k, we only show 
the exponents for two values of k, i.e. k = 0.13 
(triangles) and k = 0.16 (hexagons). The errors 
of the points are about the size of the symbols. 
In the parameter range that we have studied, the 
Liapunov exponent is roughly 0.4 — 0.5. 

4. Consequences of Irreversibility 

Due to rounding errors in the numerical in- 
tegration, the reversibility condition that is ex- 
ploited in the proof of the detailed balance con- 
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dition for molecular dynamics algorithms is vio- 
lated djUIQ. When the trajectory length r is 
increasing, we expect the strength of this viola- 
tion to grow for two reasons. The first one is due 
to the accumulation of rounding errors, which is 
expected to grow like ^/t. The second cause is 
due to the chaotic nature of the equations of mo- 
tion which amplifies the reversibility violation like 
exp(i^r) with a positive Liapunov exponent v. 

In order to quantify the effect of these viola- 
tions in a real simulation, we have run two ver- 
sions of HMC simultaneously for Wilson QCD 
with gauge group SU{2) on various lattice vol- 
umes V. One of them is with 32-bit arithmetic 
but with 64-bit arithmetic scalar products and 
summations over the lattice, which mimics the 
typical situation in simulations on a 32-bit ma- 
chine. The other one runs with complete 64-bit 
arithmetic and serves as a reference point for an 
"exact" program. 

The Metropolis step at the end of the trajectory 
in molecular dynamics algorithms depends on the 
value of AH, the difference of the Hamiltonian af- 
ter and before the integration of eq. (^) . We have 
compared the two measurements of A_ff , one ob- 
tained from the 32-bit arithmetic version and the 
other from the 64-bit version of the program. ^ 
The absolute value of the difference between the 
two values of AH is measured after each molecu- 
lar dynamics step (< \S{AH)step\ >), and at the 
end of the trajectory (< \S{AH)traj\ >), with the 
finding that < \S{AH)traj\ > to be larger than 
< \S (AH) step] >■ We also observe that the value 
of < \S{AH)step\ > is increasing linearly with 
i/i7 = while the value of < \S{AH)traj\ > 
seems to grow faster. On 8^,12"* and 16'' lat- 
tices, we find < \5{AH)traj\ > to be 0.08%, 
0.18% and 0.7% of < |AiJ| > respectively. If 
this trend is still maintained, at about a 32** size 
lattice, < \S{AH)traj\ > could reach about 5% 
of < \AH\ > for long trajectories in HMC algo- 
rithm, which we think is already dangerous. How- 
ever, it seems that the quantity < \S{AH)step\ > 
will remain reasonably small, at about 1% of AH 
even on a 32^ lattice. This suggests that running 



the Kramers equation algorithm will still be fea- 
sible on such lattices with 32-bit arithmetic. 

5. Conclusions 

We have investigated several questions related 
to the reversibility problem of the molecular dy- 
namics algorithms for the simulation of lattice 
QCD. We have determined the leading Liapunov 
exponents for both pure SU (3) gauge theory and 
full QCD for various bare parameters. We esti- 
mated that the effects of rounding errors would 
become dangerous when running the HMC algo- 
rithm on a large lattice of 0(32^) with only 32-bit 
precision. In contrast, the comparably perform- 
ing Kramers equation algorithm has substantially 
reduced rounding error effects and its use on such 
lattices is safer. 
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^ Other physical observables could also be measured and 
compared and we will report this in the near future. 



